Avalanches, branching ratios, and clustering of attractors in Random Boolean 
Networks and in the segment polarity network of Drosophila 
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We discuss basic features of emergent complexity in dynamical systems far from equilibrium by 
focusing on the network structure of their state space. We start by measuring the distributions of 
avalanche and transient times in Random Boolean Networks (RBNs) and in the Drosophila polarity 
network by exact enumeration. A transient time is the duration of the transient from a starting 
state to an attractor. An avalanche is a special transient which starts as single Boolean element 
perturbation of an attractor state. Significant differences at short times between the avalanche and 
the transient times for RBNs with small connectivity K - compared to the number of elements 
A'' - indicate that attractors tend to cluster in configuration space. In addition, one bit flip has a 
non-negligible chance to put an attractor state directly onto another attractor. This clustering is 
also present in the segment polarity gene network of Drosophila melanogaster, suggesting that this 
may be a robust feature of biological regulatory networks. We also define and measure a branching 
ratio for the state space networks and find evidence for a new time scale that diverges roughly 
linearly with N for 2 < K <^ N . Analytic arguments show that this time scale does not appear in 
the random map nor can the random map exhibit clustering of attractors. We further show that for 
K = 2 the branching ratio exhibits the largest variation with distance from the attractor compared 
to other values of K and that the avalanche durations exhibit no characteristic scale within our 
statistical resolution. Hence, we propose that the branching ratio and the avalanche duration are 
new indicators for scale-free behavior that may or may not be found simultaneously with other 
indicators of emergent complexity in extended, deterministic dynamical systems. 

PACS numbers; 05.45.-a, 89.75.-k, 89.75.Fb, 89.75.Da, 87.18.Cf 



I. INTRODUCTION 

Random Boolean Networks (RBNs) [ij are elementary 
models for signaling processes such as genetic regulation 
- where a binary state based on Boolean logic [3,[i,|3| en- 
capsulates local gene expression. The dichotomy between 
their easy construction and their emergent complex be- 
havior has motivated researchers in diverse fields includ- 
ing the neurological computational @, evolution- 
ary 0] and physical sciences to use these or related 
models as test beds for ideas about self-organization. 

The most well known fact about RBNs is that they 
exhibit three distinct phases in a statistical ensemble ob- 
tained by averaging over random realizations: chaotic, 
frozen and critical, depending on the connectivity K of 
the Boolean elements. Derrida and Pomeau used an 
annealed approximation to prove that if = 2 is a critical 
ensemble in between ordered and chaotic regimes. Ac- 
cording to their analysis, the distinct phases correspond 
to different patterns of growth for the Hamming distance 
between two nearly identical initial states. For K > 2 
the distance on average grows exponentially. For K < 2 
the distance on average decays exponentially. For K = 2 
the ensemble of random realizations is critical - the dis- 
tance is dominated by fluctuations. The existence of this 
critical phase has been used to argue that many natural 
systems, includin g lif e itself, function at a so-called "edge 
of chaos" ElilMlillll- 

However, different measures of criticality, which appear 



simultaneously in the K — 2 ensemble of RBNs, may not 
point to the same critical behavior when applied to real 
signaling networks. For instance, there may be many dis- 
tinct "edges of chaos" (according to different definitions 
of dynamical criticality) that can appear in networks that 
are not members of a statistical ensemble, but are, in- 
stead, organized by natural selection or by other forces. 
More to the point, we believe that one should not conflate 
all measures of criticality as being necessarily equivalent 
in treating complex systems far from equilibrium since 
there is no reason a priori that they should probe in- 
trinsically related dynamical fluctuations. Hence, weq 
explore several alternative probes of criticality in RBNs 
that may be more empirical or, indeed, more useful in 
revealing the forces and constraints that shape natural 
or man-made regulatory or signaling processes. 

To this end, here we develop new methods of complex 
network analysis to analyze ensembles of state space net- 
works (SSNs) for RBNs and compare their measurable 
observables to that of the segment polarity network of 
Drosophila melanogaster where possible. The as- 
pects of these SSNs that have received most attention 
so far include the probability distribution of attractor 
lengths and basins of attraction sizes [H, . In the 
context of regulatory networks, attractors are thought to 
correspond to distinct cellular states [ll| or cycles [20| . 
while the attractors of a signal transduction network cor- 
respond to steady state response(s) to the presence of 
a given signal [Sri . In Ref . [l^, [23, H^l some exact re- 
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suits were derived for RBNs with K = I. The critical 
case of K — 2 was re-exarained numerically in Ref . [2^ , 
which reported power law behavior in the distribution 
of transient times. Ref. [26| examined the probability 
of returning to the same attractor after perturbing var- 
ious number of nodes. In Ref. [13] methods of complex 
network analysis were able to distinguish the K = 2 (crit- 
ical) ensemble from the others using measures of network 
heterogeneity in the SSNs. These network measures in- 
clude node degree (a local measure) and path diversity (a 
global measure) as well as variations in the path diversity 
between different realizations in the statistical ensemble. 

Inspired by the idea of self-organized criticality p^ . 
here we start by using avalanches to probe the structure 
of the SSNs. Avalanches are the responses of a system 
in an attractor state to small perturbations (a bit flip of 
a single Boolean element). These eventually die out and 
the system returns to an attractor state, which will be 
the same or a different attractor. By exact numeration, 
we find that the distribution of avalanche times allows a 
clear distinction between the K = 2 ensemble and other 
values of K. While the avalanche and transient time dis- 
tributions converge (except for an overall normalization) 
at large times, significantly more avalanches of short du- 
rations exist compared to transients. This shows that at- 
tractors preferentially cluster in configuration space for 
RBNs for all K small compared to the number of ele- 
ments in the network, N. This feature is also found in 
the Boolean representation of the segment polarity gene 
network in Drosophila p^ . This biological network has 
a distribution of avalanche times closer to the K — 2 
ensemble than to other values of K. 

In order to clarify the differences between avalanches 
and transients at short times we define a "branching ra- 
tio" to descibe how the average number of dynamical 
states grows as a function of distance from the attrac- 
tor. The average branching ratio is given by the ratio 
of the probability distribution for transient times at suc- 
cessive times, P{Tt + l)/P{Tt). For 2 < <C TV the 
quantity crosses unity at a certain distance from the at- 
tractor where the SSN is the most "bushy" . The crossing 
time grows with system size, N , indicating a new diverg- 
ing time scale for RBNs. This scale may also separate 
the short time regime where the avalanche and transient 
times differ from the long time regime where they con- 
verge. 



A. Outline 

In Section II, we review basic facts about RBNs and 
define relevant quantities for our analysis. Section III 
contains results from numerical simulations for RBNs, 
while Section IV compares the behaviors found with that 
in the Drosophila segment polarity network. We conclude 
with a summary in Section V. 



II. RANDOM BOOLEAN NETWORKS 

An RBN consists of N Boolean (0, 1) elements where 
the value of each element evolves in discrete time accord- 
ing to a random Boolean function of K distinct input 
arguments. Each of these arguments is chosen randomly 
from the N Boolean elements of the network. For each 
set of values of the arguments, the Boolean function is 
chosen randomly to be 1 with probability (bias) p and 
with probability 1 — p. The inputs and functions as- 
signed to each element remain fixed for each realization. 
Ensemble averages are achieved here by considering many 
different realizations of the Boolean network with K and 
iV constrained to certain values and for bias, p = 0.5. 

The state of the RBN is a bit-vector that specifies 
the value of each Boolean element. Here we are exclu- 
sively concerned with the case of deterministic dynamics 
achieved by updating all N elements of a given RBN in 
parallel. This uniquely maps each state to one successor 
state, known as its "image" . Consequently, the dynam- 
ics of an RBN can be visualized as a state space network 
(SSN) [13, nil. An SSN is made by connecting each of the 
Af = 2^ dynamical states of an RBN to its image with 
a directed link. The out-degree of each node in the SSN 
is the number of its images and hence exactly one. The 
in-degree of a node, which is the number of pre-images of 
the state, ranges from to A/" in principle. In the context 
of SSNs, the "distance" between two states or nodes in 
the network is the length of the shortest path connecting 
them — if such a path exists. This directly implies that 
the distance between a state and its image is one. Al- 
ternatively, one can use the Hamming distance (number 
of different bits between two state vectors) as a metric. 
Whenever we do that, we call the set of states "configu- 
ration space" . In general, there is no simple relationship 
between these two measures of distance, i.e. between the 
configuration space and the state space. 

In the limit K = N, every Boolean element has the 
same set of inputs and its dynamics is a random function 
of the entire state of the RBN. Since the Boolean func- 
tions are chosen randomly this implies that the image for 
each state coincides with the definition of a random map. 
By construction, the associated SSN of a random map is 
a random directed graph with a Poisson distribution for 
its in-degree, with a mean of one while the out-degree of 
each node is fixed to be one. 

For finite N, any initial state eventually evolves to an 
attractor, which may be a single state or a periodic cycle. 
Depending on the specific RBN, an arbitrary positive in- 
teger of different attractors can coexist within a single 
realization of an SSN. Disconnected basins of attraction 
can occur - each of which consists of all states that evolve 
to the same attractor. Each state that does not belong to 
an attractor is called a transient state and is visited only 
once. Naturally, one can assign a quantity to each at- 
tractor cycle which characterizes the probability W that 
starting from a random node in the SSN, the system is 
found in any particular state of this attractor cycle after 



3 



an arbitrarily long time (in the steady state) . For a state 
on the ith attractor cycle, Wi is 



W, = 



J\fA, 



(1) 



where Bi is the total number of states in the basin of 
attraction and Ai is the length of the attractor cycle. 

We study the probability distribution of transient 
times and avalanche durations for ensembles of RBNs 
with fixed K and N. The transient time, Tt, for a given 
initial state is defined as the number of time steps re- 
quired to reach an attractor. States constituting attrac- 
tor cycles are assigned a transient time Tt = 0. Here, 
the distribution of transient times is obtained by consid- 
ering all possible initial states (via exact enumeration) 
over many independent realizations of the RBN. 

Avalanches, on the other hand, are created by fiipping 
a single Boolean element (0, 1 1,0) of an attractor 
state. This definition resembles avalanches in the "Game 
of Life" [3, . An avalanche continues until the sys- 
tem returns to a (possibly different) attractor. Hence, 
the avalanche durations, Tq, are the transient times for 
the collection of initial states generated by single bit flip 
perturbations to attractors. For a given RBN, the dis- 
tribution of these avalanche durations is obtained from 
all possible avalanches created by every single Boolean 
element flip of each attractor state. 

Using exact enumeration avoids many potential biases 
typically encountered when estimating properties of the 
dynamics as, for example, those related to excluding very 
long transients due to computational constraints. The 
trade-off is that we only examine relatively small systems, 
iV « 25, - although the SSNs are large 0(10^) [33]. 

To make our results comparable to other studies that 
use random sampling rather than exact enumeration ^25*1 , 
each avalanche is given a weight Wi defined in Eq. ([1]) and 
the weights are accumulated over many realizations in 
order to obtain a probability distribution for avalanche 
times in the RBN ensemble. We will discuss the effects of 
different weighting schemes in a future publication (29j . 

Finally, we define the average branching ratio as the 
ratio of the number of states (in the RBN ensemble) at 
distance T + 1 from an attractor to the number of states 
at distance T. We study the finite size properties of this 
quantity and find indications for a hitherto unknown, di- 
verging time scale for RBNs that cannot exist in the ran- 
dom map. 



III. RESULTS FOR RANDOM BOOLEAN 
NETWORKS 

Fig. [T] shows the probability distribution function 
(PDF) for avalanche durations, P{Ta), for various values 
of AT, as well as the random map. For AT = 1, P{Ta) is a 
narrow distribution. For AT > 3, a plateau in P{Ta) ap- 
pears that widens as K increases. For A' = 3 and AT = 4 
the cut-off decays exponentially while for AT > 5 and the 
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FIG. 1: (Color online). The probability distribution function 
(PDF) of avalanche times, P{Ta), for various values of K and 
the random map, all with N = 17. The K = 2 curve stands 
out as a broad distribution without an apparent cut-off in 
our measurement window. Numerical results are for 1.4 x lO'' 
realizations for A = 1, 4 x 10® for A = 2, 9 x 10^ for A = 3, 



1 x 10'' for A : 
random map. 



4, 2.6 x lO"" for A = 6 and 1.3 x lO'' for the 



random map, the cut-off decays faster than exponentially 
with T. For AT = 2, P{Ta) decays slowly with no observed 
cut-off. This is confirmed in Fig. [2l which shows the de- 
pendence of P{Ta) on N for A" = 2. While the curves 
vary with N , no characteristic time scale appears within 
our statistical resolution. This suggest that avalanche 
durations are an indicator of criticality in RBNs. 

Fig. [2] also compares P{Ta) and P{Tt) for AT = 2. 
While P{Ta) and P{Tt) clearly differ for small arguments, 
the two distributions become statistically indistinguish- 
able, up to an overall normalization, for larger arguments. 
This is actually true for all values of AT studied here — 
see Fig. |3]for the case AT = 6. For larger K , the curves 
approach each other more quickly and for K = N the 
two PDFs are statistically and theoretically identical. As 
mentioned previously, the K = N case corresponds to 
the random map, where no correlation exists between the 
Boolean values of the elements in a state and its image. 
As a result, flipping a single bit on an attractor state can 
put the system into a state anywhere in the entire state 
space, and the distribution of avalanche times, P{Ta), 
for the random map is identical to the distribution of the 
transient times, P{Tt). 

Ref. [25] claimed that the PDF for transient times for 
AT = 2 RBNs has a power-law tail for large N. Based 
on our results shown in Fig. 2 this also suggests that 
the distribution of avalanche times could have a power- 
law tail for large N reminiscent of self-organized critical 
systems [l2|. Note, however, that the purported "expo- 
nent" of the power-law decay in Ref. [25| varies system- 
atically with system size (as shown in Fig. [5]) and that 
the behavior in the limit of infinite system size remains 
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FIG. 2: (Color online). Avalanche and transient times for 
K ^ 2 RBNs with diflFerent N. PDFs for avalanche times 
are indicated by solid lines. Dashed lines correspond to the 
PDFs, P{Tt), for transient times. These numerical results for 
avalanche (transient) times are based on the following number 
of realizations: 7.5 x 10*^ (9 x 10^) for N ^7, lO'' (9 x lO'') 
for iV = 11, 4 X 10*^ (2.8 x 10®) for N = 17, and 2.5 x lO'^ 
(2.6 x 10*) for iV = 23. Note the differences at short times 
while the long time behavior is statistically identical up to an 
overall normalization. 



FIG. 3: (Color online). Avalanche and transient times for 
K = Q RBNs. The PDFs for the avalanche (transient) times 
are indicated by solid (dashed) lines. The results for avalanche 
(transient) times are based on the following number of real- 
izations: lO'^ (8 X 10") for Af = 7 - 15, 2 X 10* (2 x 10*) for 
N = 19, and 6 x 10^ (250) for N = 23. Note that as for 
K = 2, also for K — 6 there are differences at short times, 
while the long time behavior is statistically the same up to 
an overall normalization. 



unclear. One possibility is that (as discussed in detail 
later) the position of maximum in the transient time dis- 
tribution increases roughly with N while the largest pos- 
sible time is bounded by the number of states N — 2^ . 
This would give an asymptotically uniform distribution 
at large times for large N . 

The aforementioned differences between the PDFs for 
avalanche times and transient times are characterized by 
P{Ta) > P{Tt) for small arguments for all studied val- 
ues of N with K < N . Both the transient time and 
avalanche duration are the time it takes for the system 
to reach an attractor from an initial state. While for the 
transient times the initial state is chosen randomly, for 
avalanche times the initial state is chosen by making a 
perturbation to an attractor state. We conclude that a 
state generated by a single flip perturbation to an attrac- 
tor state tends to be closer to an attractor state in the 
SSN than a randomly chosen state. It also has a statis- 
tically significant chance to be, itself, an attractor state 
on a different attractor. 

Fig. 0] shows the cumulative distribution (CDF) of 
avalanche durations for all perturbations, for perturba- 
tions that remain in the same basin of attraction and for 
those that lead to a different basin of attraction and com- 
pares them with the CDF of transient times, for K = 2. 
A single element perturbation leads to a nearest neighbor 
in configuration space. Fig. [3] shows that these neighbors 
in configuration space are more likely to be close to a 
(possibly different) attractor in the SSN than randomly 
selected states. Fig. U also indicates that there is a non- 
negligible probability that avalanches that go to different 



attractors have duration zero. This is indicated by the 
fact that the curve c does not go through the point (0,0). 
It means that one bit flip is sufficient to put an attractor 
state directly onto a different attractor. This suggests in 
particular the existence of regions in configuration space, 
where the attractors are more likely to live: Different at- 
tractors in the SSN are often close in configuration space. 

Qualitatively similar results hold for all K ^ N. How- 
ever, independent of N, the fraction of avalanches that re- 
turn to the same attractor monotonically decreases as K 
increases, starting at 95% for K = 1 and asymptotically 
converging to 2/3 for the random map for large N [26| . 
Moreover, the PDF for the duration of avalanches that 
lead to different attractors as well as the corresponding 
PDF for avalanches leading to the same attractors ap- 
proach the PDF for transient times as K approaches N. 

To see why the attractors are clustered in configuration 
space, we consider the notion of relevant and irrelevant 
components (for a thorough discussion see e.g. [sol. IsH). 
The elements of an RBN can be divided into relevant 
and irrelevant elements. The relevant elements are clus- 
tered into what are called relevant components. These 
components determine the number and lengths of attrac- 
tors [sol . [sH . The irrelevant elements are those whose 
information is lost in the dynamics. As an example, con- 
sider the simple network in Fig. [S) The network is made 
of the relevant component made of elements A and B and 
the irrelevant component element C. A and B form an 
information conserving loop by copying each other, while 
C is forced to 1 regardless of the state of its inputs. More- 
over, C does not influence the state of A and B. Consider 
the attractor state [a{A), a{B), a-{C)] — [0,0,1], where 
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FIG. 4: (Color online). A comparison of the cumulative distri- 
bution functions (CDFs) in if = 2 RBNs for: (a) the duration 
of all avalanches (b) the duration of avalanches that return to 
the same attractor (c) the duration of avalanches that lead to 
a different attractor (d) the transient times. The inset shows 
a linear magnification of the upper left section of the original 
plot. Curves are based on the following numbers of realiza- 
tions of RBNs with N = 17: (a) on 4 x 10*^, (b) and (c) on 
7 X lO'^ and (d) on 2.8 x 10^ 



a{X) is the value of the element X. If we flip element C 
the system will return to the same attractor [0, 0, 1] once 
the information injected into the irrelevant component is 
lost. Alternately we could flip element A which is part of 
the relevant component. This would lead to state [1, 0, 1], 
which is a state on a period two attractor cycle composed 
of [1,0,1] and [0,1,1]. This perturbation leads directly 
to a different attractor. In general we can replace the 
element C by a set of elements {Ci} which are neither 
influenced by nor influence A and B. This allows us to 
the see the clustering of attractor states; the attractor 
states: [0, 0, a({Ci})], [0, 1, a({CJ)], [1, 0, a({Q})] and 
[1, 1, tT({Ci})] are clustered into four neighboring sites in 
configuration space. 

In general flipping the state of an irrelevant element 
on an attractor will be forgotten and the dynamics even- 
tually leads to the same attractor. However, by flipping 
the state of a relevant element, the system can reach a 
different basin of attraction. 

In the chaotic phase, K > 2, the distribution of tran- 
sient times can be estimated from the joint probability 
distribution for transient times and attractor lengths de- 
rived in [3^ ] using an annealed approximation, 



^ 2" -Tt 



A=l 



exp 



1 f {Tt+A) 



(Tt + l)/(e°"/2) 



(2) 



for large N and T^. In general, a depends on K. It was 
shown in Ref. [32] that a vanishes for K = 2, and ap- 
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FIG. 5: The network structure for an illustrative example of 
an RBN with K — 2. Elements A and B copy each other at 
each time step while element C is forced to 1. 




FIG. 6: (Color online). The PDF P{Tt) for the duration 
of transients in random maps (used to simulate the SSNs of 
K = N RBNs) for various A'' collapsed using Eq. ((2]). The 
inset shows the un-scaled PDFs. The numerical results for 
= 7 — 17 are based on 1.4 x 10^ realizations of the random 
map, A = 19 on 3.4 x 10^ and A^ = 21 - 23 on 500. 



proaches the value log 2 0.69 for the random map. Bas- 
toUa and Parisi also pointed out that, while the analytical 
expressions for the transient time distribution may not be 
good approximations for small K, the time scale e"-^!"^ 
derived using the annealed approximation (see Eq. ^) 
is in good agreement with their simulations for if = 3. 

This last result is conflrmed by our numerical simula- 
tions for various values of K. Fig. [5] shows that P{Tt) for 
the random map and various N indeed exhibits finite size 
scaling with characteristic time 2^/^. For if = 6 we also 
find a reasonable data collapse, shown in Fig. [7| using the 
value a — 0.58 obtained in [S^- Yet, the quality of the 
data collapse continues to deteriorate as K decreases to 
3. 
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FIG. 7: (Color online). The PDF P{Tt) for the duration of 
transient times in A' = 6 RBNs collapsed using Eq. ([2}. The 
inset shows the un-scaled PDFs. The numerical results for 
= 7 — 15 are based on 8 x 10'' realizations of the random 



map, iV = 17on 2.5 x 10^ 
on 250. 
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increase becomes more prominent for smaller K. 

The number of states in the ensemble that are at a dis- 
tance Tt from an attractor is 2^P{Tt). Considering the 
arboreal structure of the SSN, this number corresponds 
to the number of states in the Jt-th shell. An increase 
of P{Tt) for small Tt then means that the shells become 
more populated with distance from an attractor. The 
monotonic decrease of Eq. ^ indicates that the Tt = 
(attractor) shell should be the most populated, which we 
do observe for the random map. The maxima observed 
in Fig. [5] for RBNs with K < N indicate the existence of 
a new characteristic distance for RBNs with K < N that 
does not appear in the random map. 
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FIG. 9: (Color online). The mean branching ratio, P(Tt + 
l)/P{Tt), for the arboreal structure of the SSNs for RBNs for 
various values of K and the random map with A'^ = 17. The 
inset shows a magnification of the same curves. The mean 
branching ratio of 1 — which is the mean degree of the SSN 
— is shown for comparison. The results are for 2.8 x 10® 
realizations for if = 2, 8 x 10^ for = 3, 9 x lO'^ for iC = 4, 
2.5 x lO'^ for K = 6, and 1.3 x lO'^ for the random map. 



FIG. 8: (Color online). The PDF, P(Tt), for transient lengths 
in the SSNs of RBNs for various values of K and the random 
map, A'^ — 17. Each (excluding K = 2) is accompanied by a 
corresponding theoretical curve based on Eq. ([2]) shown with 
dashed line. Numerical results are for 2.8 x lO*" realizations 
for it" = 2, 8 X 10^ for it = 3, 9 x 10^ for K = A, 2.5 x lO'^ for 
K — 6, and 1 x lO"* for the random map. Based on [3^ . we 
used the values a = 0.20, 0.38, 0.58 and log 2 for if = 3, 4, 6 
and the random map respectively. 



Fig. [8] shows the distribution of transient times along 
with the theoretical predictions given by Eq. ^ . For the 
random map, the prediction is exact and the agreement 
between data and theory is excellent. For K = 6 the the- 
oretical predictions do not match the data for small Tt- 
The theoretical curve in Eq.(I2]) decreases monotonically 
with T for all K > 2. However, the actual distribution 
increases initially with Tt for all i^T < considered. This 



The ratio of the sizes of consecutive shells are shown 
in Fig. [HI A state in the (Tt -I- l)-th shell can be seen as 
branching from its image in the Tt-th shell. If each basin 
of the SSN is viewed as a tree with the attractor at its 
root. Fig. [3] gives a mean branching ratio from the Tt-th 
shell to the (Tt + l)-th shell for the whole ensemble. Note 
that this is not the mean of the branching ratios for each 
attractor basin. For K < N, Fig. [9] shows that for small 
Tt the shells first grow at a rate that slows down with dis- 
tance from the attractors. The growth eventually stops 
and the shells become smaller and smaller for increasing 
Tt when it is larger than a certain value that depends on 
K (and on N). This turn around happens when the ratio 
drops below unity. For the random map, Eq. ^ implies 
that shells should become smaller for all Tt, albeit very 
slowly. This can be seen in Fig. [9] in which the random 
map curve almost coincides with unity. The K — 6 curve 
crosses unity and quickly merges with the random map 
curve, while the if = 4, 3 and 2 curves make a discernible 
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dip below the P{Tt + l)/P{Tt) = 1 line, the extent of the 
dip depending on K. Indeed it appears that the behav- 
ior of the branching ratio for K — 2 RBNs shows more 
variation than other values of K. 




FIG. 10: (Color online). The value of Tt where the mean 
branching ratio (Fig.[9)l is minimal and the value where P{Tt + 
l)/P{Tt) = 1 for various sizes of RBNs with K = 2. 

Both the value of Tt where the branching ratio is 
minimal and the value where P{Tt + l)/P{Tt) = 1 in- 
crease with increasing N as shown for K — 2 m Fig. [TOl 
The roughly linear increase also appears to hold for 
2 < K < N. If the growth that we find in Fig. [TO] can 
be extrapolated to large system sizes, our results indicate 
the existence of a new diverging time scale for RBNs with 
2 < K <^ N that is unrelated to the cut-off for large ar- 
guments given by Eq. 

IV. RESULTS FOR THE SEGMENT POLARITY 
NETWORK OF DROSOPHILA MELANOGASTER 

To compare our results for RBNs with a biological sig- 
naling network, we consider the segment polarity network 
of fruit flies [l^ . This Boolean model of the gene and pro- 
tein interactions involved in embryonic pattern formation 
in the fruit fly Drosophila melanogaster is a particularly 
well-documented and successful application. In particu- 
lar, it reproduces many experimentally observed features, 
including knock out results, indicating that knowledge of 
the kinetic details are not necessary. 

The model presented in [Tg'I considers gene expression 
patterns in four adjacent cells that form a Drosophila 
parasegment. Each cell is modelled by a network of flf- 
teen Boolean elements. Each element represents either 
an mRNA species or a protein species in the cell. The 
state of each element is one if the corresponding species 
is expressed, otherwise it is zero. Each element in a cell 
can interact with elements within the same cell as well 
as elements in the neighboring cells. All interactions are 
modelled by Boolean functions. Of the 15 species, the 



state of protein "sloppy-paired" (SLP) in each of the 4 
cells is fixed to its biologically relevant value. The result- 
ing state space contains 6 fixed point attractors. 

Here we study the dynamics of only one cell, for two 
reasons. First, the complete network, which consists of 
60 elements, is computationally untractable using exact 
enumeration. Second, a single cell represents the build- 
ing block of the entire organism. Understanding the dy- 
namics of a single cell is important to understand the 
whole network. Robust dynamics of the individual cells 
that form the segment has been argued to underlie the 
observed robustness of expression patterns of the whole 
segment [ssj .A differential equation model for individual 
cells based on the original work in [s^l was used to show 
this relation between the dynamics of the individual cells 
and the whole segment. In the same vein, here we study 
the state space of the Boolean network for an individual 
cell. As a result, our individual cell network cannot be 
directly compared with the parasegment network studied 
in 0. 

In the Boolean network model of [l^ , cells interact via 
three elements, the hedgehog mRNA (hh), the hedgehog 
protein (HH), and the wingless protein (WG). To reduce 
the model to one individual cell, we replace the values 
of these three elements in the cell's neighbourhood by 
three new variables, hhN, HHN and WON. Each of these 
three elements evolves by copying its state. Furthermore 
SLP also copies its own state. Therefore, the values of 
the states of these four elements are fixed by the initial 
conditions and do not evolve dynamically. Note that this 
ensures the existence of at least 16 different attractors, 
each corresponds to a unique combination of the input 
elements. 

The polarity network studied here contains 15 + 3 = 18 
elements in total. The connectivity, K, of the nodes in 
the polarity network is distributed as follows: 9 nodes 
have K — 1, A have K = 2, 4 have K = 3 and only one 
node has K = A. The mean connectivity in the network 
is thus 1.83. The Boolean rules (frmctions) of the nodes 
in the network are given in Table HI The average bias, p, 
in the network is 0.42. The internal homogeneity, which 
measures the average percentage of the most abundant 
outcome (zero or one) in each Boolean function, is 0.625. 

We find that the individual cell state space has 21 dif- 
ferent attractors, each of which is a fixed point. The 
states of the 21 attractors are shown in Fig. [TT] Each of 
21 columns in the figure represents an attractor. Each 
of 18 rows represent the value of the Boolean state of a 
particular element across different attractors. The black 
squares depict ones while white squares depict zeros. 

Fig. [TT] also illustrates the clustering of attractors in 
configuration space. Attractors in the figure are arranged 
in clusters where neighboring attractors in each cluster 
are one bit fiip away. Note that the bit flips that im- 
mediately lead to a different attractor are all flips of the 
relevant component. In this context, there are at least 
4 relevant components in the network SLP, hhN, 

HHN, and WGN all copy their own state at each time 
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Node Boolean function (input/output relation) 

^^g II (CIA AND SLP AND NOT CIR) 

OR (wg AND (CIA OR SLP) AND NOT CIR) 

WG wg 




WGN AND NOT SLP 



EN 


en 


hh 


EN AND NOT CIR 


HH 


hh 


ptc 


CIA AND NOT EN AND NOT CIR 


PTC 


ptc OR (PTC AND NOT HHN) 


PH 


PTC AND HHN 


SMO 


NOT PTC OR HHN 


ci 


NOT EN 


CI 


ci 


CIA 


CI AND (SMO OR hhN) 


CIR 


CI AND NOT SMO AND NOT hhN 


SLP 


SLP (input) 


hhN 


hhN (input) 


HHN 


HHN (input) 


WGN 


WGN (input) 



TABLE I: The Boolean function of each node in the 
Drosophila segment polarity network. Upper case names are 
reserved for proteins and the lower ceise names are for mRNA. 
The hhN, HHN and WGN nodes refer to the state of hh, HH 
and WG in the neighboring cell. 



step, therefore each is a relevant component. 

We note that the attractors in Fig. [11] capture the cell 
states observed in Drosophila segments. For example, 
none of the cells in the Drosophila segment co-express 
genes wg and en [3^1 ■ Lack of wg-en co-expression is also 
seen in the attractors in Fig. [TTJ Our results suggest that 
the observed gene expression phenotype of the segment 
polarity network stems from the constrained dynamics 
in a single cell, rather than the complex interactions be- 
tween the cells. The same has been argued previously 
in [3^ using a differential equation model for a single 
cell. The biological significance of these results and their 
impact on gene expression robustness will be expounded 
in a future publication [36j . 

Avalanches and transients for the segment polarity net- 
work are defined in the same way as for RBNs — see 
section |lTl Fig. [T^] shows the distributions of avalanche 
and transient times for the segment polarity network in 
Table H] As for RBNs (see, for example. Fig. H]) there 
are more short avalanches than short transients indicat- 
ing that attractors cluster in configuration space. This 
holds for the specific avalanches that return to the same 
attractor as well as those that go to a different attractor. 
Interestingly, in the latter case there is a significant prob- 
ability that an avalanche will immediately reach another 
attractor — similar to the case of RBNs with K ^ N. 

The qualitative similarities between RBNs and the 



FIG. 11: The attractor states of the Drosophila polarity net- 
work. Each column represent an attractor, the dark and white 
square represent 1 and respectively. The order of the attrac- 
tors was chosen to show attractor clustering in configuration 
space. 




FIG. 12: (Color online). Comparison of CDFs from the seg- 
ment polarity network of Drosophila (N — 18) for: (a) the 
duration of all avalanches (b) the duration of avalanches that 
return to the same attractor (c) the duration of avalanches 
that lead to a different attractor (d) the transient times. The 
inset shows the PDFs for the same data. 



Drosophila polarity network are also demonstrated in 
Fig. [itl where the distributions for RBNs with K =1,2,3 
are shown for comparison. The large scatter, due to the 
fact that we are examining only one network rather than 
an ensemble, prevents any precise quantitative compari- 
son between the segment polarity network and RBNs. It 
is interesting to note, though, that the distributions of 
avalanche durations for the polarity network have a ten- 
dency to be more similar to RBNs with K = 2 at small 
values of T than any other K. 
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FIG. 13: (Color online). Comparison of CDFs from the seg- 
ment polarity network of Drosophila {N = 18) with those for 
RBNs with A'^ = 18 and K = 1,2 and 3. For each case we 
show the results for (a) the duration of all avalanches and (b) 
the transient times. Note that the data for Drosophila cor- 
respond to a single Boolean network while the data for the 
RBNs correspond to ensemble averages. 



V. CONCLUSIONS 

We have studied avalanches on Random Boolean Net- 
works and on the segment polarity network of Drosophila 
by performing bit-flip perturbation to attractor states 
and waiting for the dynamics to relax back to an attrac- 
tor. In both cases (assuming K <^ N for RBNs), there 
are more avalanches of short duration than what would 
be expected based on the distribution of transients times. 
This indicates that attractors tend to cluster in configu- 
ration space. 

For the random map, which corresponds to RBNs with 
K = N, the distribution of avalanche durations, P{Ta), 
is identical to the distribution of transient times, P(Tt). 
As analytically confirmed, the distribution tends to form 
a plateau with a sharp finite-size cut-off. For 2 < K < N , 
deviations from this behavior occur for small arguments 
such that the plateau eventually disappears as K 2. 
In particular, the deviations for small arguments are dif- 
ferent for P{Ta) and P(Tf). The critical case K ^ 2 
exhibits different behavior. In this case, both distribu- 



tions are broad and neither the plateau nor the cut-off are 
observed. Their scale-free appearance suggests that the 
distribution of avalanche durations as well as the distri- 
bution of transient times can be used as possibly distinct 
indicators for the criticality in discrete, deterministic dy- 
namical systems. Indeed we find that the avalanche dura- 
tions for the Drosophila segment polarity network follows 
more closely the K = 2 behavior (compared to other val- 
ues of K) while no clear statement can be made about 
transient times for the Drosophila network we analyzed. 

The similarity of P{Ta) and P{Tt) for large arguments 
and K > 2 indicates that initial avalanche states that 
are far away from an attractor are independent of the 
states on the attractors (this is true for all arguments in 
the case of the random map). The arboreal structure of 
the SSN, the mean branching of which was presented in 
Fig[9l also seems to become more similar to that of the 
random map for large K. 

However, the branching ratio is a non-monotonic func- 
tion of Tf, crossing unity and then passing through a 
minimum before it starts to resemble the random-map 
behavior for larger Tf . Along with our results for the dif- 
ferences between P{Ta) and P{Tt) for small arguments, 
this indicates that while states lying closer to the attrac- 
tors in the SSN are correlated with the attractor states, 
this correlation decays as one moves away from the at- 
tractors on the SSN. The time scale at which this cor- 
relation in the arboreal structure of the SSN decays, or 
in which the growth rate of the tree structure in state 
space changes from increasing to decreasing may be a 
new characteristic scale for RBNs. It appears to diverge 
roughly linearly with the system size N . This time scale 
does not appear in the random map and is not related 
to that previously studied in [32], which grows exponen- 
tially with the system size N . Beyond that scale, the SSN 
of an RBN in the non-frozen regime may correspond to 
a random map, while the regions around attractors are 
distorted by their presence. 
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